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Recent experimental studies focusing on the morphological properties of surfaces 
eroded by ion-bombardment report the observation of sclf-affine fractal surfaces, 
while others provide evidence about the development of a periodic ripple struc- 
ture. To explain these discrepancies we derive a stochastic growth equation that 
describes the evolution of surfaces eroded by ion bombardment. The coefficients 
appearing in the equation can be calculated explicitly in terms of the physical 
parameters characterizing the sputtering process. Exploring the connection be- 
tween the ion-sputtering problem and the Kardar-Parisi-Zhang and Kuramoto- 
Sivashinsky equations, we find that morphological transitions may take place when 
experimental parameters, such as the angle of incidence of the incoming ions or 
their average penetration depth, are varied. Furthermore, the discussed methods 
allow us to calculate analytically the ion-induced surface diffusion coefficient, that 
can be compared with experiments. Finally, we use numerical simulations of a one 
dimensional sputtering model to investigate certain aspects of the ripple formation 
and roughening. 



1 Introduction 

In the last decade we have witnessed the development of an array of theoretical 
tools, ideas and techniques intended to descxibfi^nd characterize the growth 
and roughening of nonequilibrium surfaces El'ErHcl. Initiated by advances in 
understanding the statistical mechanics of various nonequilibrium systems, it 
has been observed that for most surfaces in nature the roughness follows simple 
scaling laws. These surfaces are self-affine fractals, being characterized by 
the roughness or self-affine exponent a. One of the main advantage of this 
description is that various growth processes can be classified into universality 
classes that share the same scaling exponents. On the practical side this means 
that the scaling exponents characterizing roughness do not vary continuously, 
but are defined by the universality class to which they belong. 

One particularly important thin film processing technique is ion beam 
sputtering n un. Sputtering is the removal of material from the surface of solids 
through the impact of energetic particles. It is a widespread technique, used 
in a large number of applications, with a remarkable level of sophistication. It 
is a basic tool in surface analysis, depth profiling, sputter cleaning, microma- 
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chining, and sputter deposition. 

Motivated by the advances in understanding growth, and by the need of 
having a detailed knowledge on the morphology of the sputter eroded surfaces, 
recently a number of experimental studies have investigated the morphological 
properties of surfaces eroded by ion bombardment. Briefly, the experimen- 
tal results can be classified in two main classes. There exists ample evidence 
ab out- Ihe. pdmielapment of a periodic ripple structure in sputter etched sm- 
faces ffl-HllilO'ES. However, a number of recent investigations have provided 
rather detailed and convincing experimental evidence, that under certain ex- 
perimental conditions ion-eroded surfaces are rough and self-affine, and the 
roughness follows the predictions of various scaling theories liSlla^ta. Moreover, 
these investigations did not find evidence of ripple formation on the surface! 

The discrepancy between the results of the mentioned investigations moti- 
vated us to have a .second look at the mechanisms shaping the morphology of 
ion eroded surfaceslI3. In this paper we investigate the large scale properties of 
ion-sputtered surfaces aiming to understand in an unified framework the var- 
ious dynamic and scaling behaviors of the experimentally observed surfaces. 
For this we derive a stochastic nonlinear equation that describes the time evo- 
lution of the surface height. The coefficients appearing in the equation are 
functions of the physical parameters characterizing the sputtering process. We 
find that transitions may take place between various surface morphologies as 
the experimental parameters (e.g. angle of incidence, penetration depth) are 
varied. Namely, at short length-scales the equation describes the development 
of a periodic ripple structure, while at larger length-scales the surface mor- 
phology may be either logarithmically (a = 0) or algebraically (a > 0) rough. 
Furthermore, we calculate analytically the ion-induced diffusion constant, D 1 , 
and its dependence on the ion energy, flux, angle of incidence, and penetration 
depth. We find that there exists a parameter range when ion bombardment 
generates a negative surface diffusion constant, leading to morphological insta- 
bilities along the surface, affecting the surface roughness and the ripple struc- 
ture. The effect of ion-induced diffusion on the morphology of ion-sputtered 
surfaces is summarized in a morphological phase diagram, allowing for direct 
experimental verification of our predictions. Finally, we use numerical simula- 
tions of a one dimensional sputtering model to investigate certain aspects of 
the ripple formation and roughening. 

2 Scaling theory 

A common feature of most non-equilibrium rough interfaces 0ffl observed ex- 
perimentally or in discrete models is that their roughening follows simple scal- 
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ing laws. The associated scaling exponents can be obtained using numerical 
simulations or stochastic evolution equations. The morphology and dynam- 
ics of a two-dimensional rough surface can be characterized with the interface 
width, defined by the rms fluctuation in the height variable h(x, y, t), 



W(L,t)= ±( J2 [h(x,y,t)-h]*), (1) 

where L is the linear size of the sample, the brackets (...) denote ensemble 
average, and the mean height of the surface, h, is defined by 

^=j2 E h &y,t). (2) 

x,y— 1,L 

For times t ^> t x ~ L z , the surface width behaves as 

W{L,t)~L a , (3) 



where a is the roughness exponent and z is the dynamic exponent. Regarding 
the early dynamics of the roughening process, the total width increases as 
W(L,t) ~ t@ , where (3 is the— growth exponent. The dynamic exponent is 
related to a and (3 as z — a//3cM. 

To understand the roughening process, we need to develop methods to 
predict the value of the scaling exponents a and (3. A breakthrough in this 
direction was the introduction of the Kardar-Parisi-Zhang (KPZ) equationliia 

dh 

— = v V 2 h + \{Vh) 2 + r 1 (x,y,t) [KPZ] . (4) 
at 

The first term on the rhs describes the relaxation of the interface due to the 
surface tension v and the second is a generic nonlinear term incorporating 
lateral growth. The noise, r](x,y,t), reflects the random fluctuations in the 
growth process and is an uncorrelated random number that has zero configu- 
rational average. For one dimensional interfaces the scaling exponents of the 
KPZ equation are known exactly, as a = 1/2, — 1/3, and z = 3/2. However, 
for higher dimensions they are known only from numerical simulations. For 
the physically most relevant two dimensional interface we have a ~ 0.38 and 
/3-0.18E3. 

If A = in (^), the remaining equation describes the equilibrium fluctua- 
tions of an interface which tries to minimize its area. This equation, introduced 
and studied in the context of interface roughening by Edwards and Wilkinson 
(EW) Ell, can be solved exactly due to its linear character, giving the scaling 
exponents a = (2 — d)/2 and = (2 — d)/4. For two dimensional interfaces we 
have a — [3 = 0, leading to a logarithmic roughening of the interface. 
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3 Experimental results 

The morphology of surfaces bombarded by energetic ions has long fascinated 
the experimental community. Lately, with the development of high resolution 
observation techniques, this question is living a new life. 

We shall focus here on two dominant morphologies, ripple formation and 
kinetic roughening, since these are observed in the sputtering of impurity free, 
amorphous materials. Impurities that bind strongly to the surface (being thus 
difficult to sputter) may [induce other dominant morphological features, such 
as cones or abrupt walls u. These will not be considered in this paper. Also, 
we limit ourselves to sputtering by ion bombardment, in which the ions have 
parallel trajectories and the same velocity. Thus we will not consider plasma 
etching (where the ions have a broad energy distribution and random angles of 
incidence) or chemical sputtering, where the yield is influenced by the chemical 
reactions taking place on the surface. 

3.1 Ripple formation 

Ripple formation on ion-sputtered surfaces have been_observed by many groups 
in various systems and ion beams (for a review see@). Here we discuss a few 
recent investigations that characterized in great detail the observed morpholo- 
gies. 

Evidence for the ripple structure on the surfaces o,f Si02 and Ge has been 
provided in a series of studies by Chason et al. We shall discuss 

here the results obtained on SiC>2 t30. A low energy ion beam (Xe, H or 
He), with energies < 1 keV is directed towards a SiC>2 sample with an angle of 
incidence of 55° from normal. The typical incoming flux is 10 13 cm _2 s _1 . The 
interfaces are analyzed using in situ energy dispersive x-ray reflectivity and ex 
situ atomic force microscopy (AFM). Bombarding the surface with 1 keV Xe 
ions, one finds that the interface roughness, determined from X-ray diffraction, 
increases linearly with the fluence (the fluence is the number of incoming atoms 
per surface area, and plays the role of time in these measurements). Thus 
[3 = 1, too large a value to be interpretable by continuum theories. Such a large 
value of [3 indicates the existence of an instability in the system. Physically, 
the instability is balanced by surface diffusion, leading to the appearance of the 
ripple structure whose wavelength increases with temperature. Such a ripple 
structure can be seen if one inspects the AFM pictures of the interface. A 
similar ripple structure has been observed for Ge surfaces bombarded by Xe 
atoms □. 

Another series of experiments on ripple formation were reported by Ma- 
cLaren et al. lid. They studied InP and GaAs bombardment with 5 keV Ar + , 
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17 keV Cs + and 5.5 keV Oj" beams in a temperature range of —50 to 200 
°C. Their study revealed in detail the temperature dependence of the ripple 
wavelength. For example for GaAs bombarded by Cs + ions the ripple spacing 
increased from zero 0.89 jxm to 2.0 /im, as the temperature increased from °C 
to 100 °C. Probably the most interesting finding of their study was that when 
lowering the temperature, the ripple spacing (wavelength) did not go contin- 
uously to zero, as one would expect, since the diffusion constant decreases 
exponentially with the inverse temperature, but rather at around 20 °C it 
stabilized at an approximately constant value. MacLaren et al. interpreted 
this as the emergence of a radiation enhanced diffusion, that gives a constant 
(temperature independent) contribution to the diffusion constant. We shall 
return to ion enhanced diffusion in Sects. 6. and 7. 

3.2 Roughening 

For graphite bombarded with 5 keV Ar ions, Eklund et al. reported a ~ 
0.2 — 0.4, and z ~ 1.6 — 1.8, values consistent with the_pcedictions of the 
Kardar-Parisi-Zhang (KPZ) equation in 2+1 dimensions li§E3. No trace was 
found of a periodic ripple structure. In these experiments pyrolytic graphite 
was bombarded by 5 keV Ar ions, which arrived with an angle of incidence 
of 60°. The experiments were carried out for two flux values, 6.9 x 10 13 and 
3.5 x 10 14 ions/cm 2 , and the total fluences obtained were 10 16 , 10 17 and 10 18 . 
The etched graphite was examined using STM. Large scale features develop 
with continuous bombardment, the interface becoming highly correlated and 
rough. 

A somewhat larger roughness exponent has been measured for samples 
of iron bombarded with 5 keV Ar arriving with angle of incidence of 25°. 
The interface morphology was observed using STM, and the height=height 
correlation function results in a roug hness exponent a = 0.53 ± 0.02 El The 
mechanism leading to such a roughness exponent is not yet understood in 
terms of the continuum theories, since for two dimensions the growth equations 
predict 0.38, 2/3 and 1, all values far from the observed value. 

Finally Si(lll) sputtered by 0.5 keV Ar + ions has also been observed to 
roughen, in this case following_an anomalous dynamic scaling form w(t,£) ~ 
ln(t)£ 2a , with a ~ 1.15 ± 0.08ll3, where w(t, £) is the surface width for a small 
window of lateral size L 

4 Continuum theory 
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4-1 Sigmund's theory of sputtering 

In order to calculate the sputtering yield, and predict the surface morphology, 
we first need to understand the mechanism of sputtering, resulting from the 
interaction of the incident ion and the surface layer. 




Figure 1: Reference frames for the computation of the erosion velocity at point O. Inset: 
Following a straight trajectory (solid line) the ion penetrates an average distance a inside 
the solid (dotted line) after which it completely spreads out its kinetic energy. The dotted 
curves are equal energy contours. Energy released at point P contributes to erosion at O. 

A qualitative picture is as followS (see Fig. Q). The incoming ions pene- 
trate the surface and transfer their kinetic energy to the atoms of the substrate 
by colliding with the substrate atoms, or through other processes such as elec- 
tronic excitations. Atoms that recoil with sufficient energy undergo secondary 
collisions, thereby generating another generation of recoiling atoms. A vast 
majority of atoms will not gain enough energy to leave their lattice position 
permanently. However, some are permanently removed from their sites, lo- 
cally making the substrate amorphous. The atoms that are near the surface 
and gain enough energy to break their bonds and leave the surface will be sput- 
tered. The scattering events that might lead to sputtering take place within 
a certain layer of average depth a. Usually the number of sputtered atoms is 
orders of magnitudes smaller than the total number of atoms participating in 
the collision cascade. 

A rather successful theory of the above processes was introduced by Sig- 
mund to describe the experimentally observed sputtering yields E3. His treat- 
ment considers the energy transfer from the incoming ion to the atoms of an 
isotropic solid by writing down a Boltzmann transport equation for the atoms. 
Expanding this equation in form of Legendre polynomials, he obtains a solu- 
tion using the method of moments. One of the most important result of his 
analysis is that for low energies the damage and energy distribution generated 
by the incoming ion follows a Gaussian. Thus here, followingE3l23, we consider 
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that the average energy deposited at point O due to the ion arriving at P 
follows the Gaussian distribution 



^ ) - (27r) 3/ 2ff „2 ex P (-^5 - -^?-) • (5) 

In (^|) z' is the distance measured along the ion trajectory, and x' , y' are 
measured in the plane perpendicular to it (see Fig. for simplicity in the figure 
x' has been set to 0); e denotes the total energy carried by the ion and a and 
fj, are the widths of the distribution in directions parallel and perpendicular 
to the incoming beam, respectively. However, the sample is subject to an 
uniform flux J of bombarding ions. A large number of ions penetrate the solid 
at different points simultaneously and the velocity of erosion at O depends on 
the total power So contributed by all the ions deposited within the range of 
the distribution (0). If we ignore shadowing effects among neighboring points, 
as well as further redeposition of the eroded material, the normal velocity of 
erosion at O is given by 

v=p / dr $(r) E(r), (6) 
Jn 

where the integral is taken over the region TZ of all the points at which the 
deposited energy contributes to So, 3>(r) is a local correction to the uniform 
flux J and p is a proportionality constant between power deposition and rate 
of erosion. In the following we revieaLiJifi .basic steps in the calculation of v; 
further details can be found in Refs. EIrEJcJ. 

4-2 Continuum equation for the surface height 

In this section we derive an equation of motion for the surface height from the 
physical model of ion-sputter erosion discussed in the previous section. Since 
we are mainly interested in the physically relevant case of a two dimensional 
substrate and the one dimensional case to linear order is very clearly explained 
in the work by Bradley and Harper c3, we refer the reader to that reference, 
and focus here on the more general 2d case. 

In the following we summarize the steps in the derivation of the equation 
of motion. 

(i) First we calculate the normal component of the velocity of erosion vo 
at a generic point O of the interface. This calculation is most easily performed 
in a local frame of reference (X, Y, Z) defined as follows: the Z axis is identified 
with the normal direction to the average surface orientation at O. Moreover Z 
forms a plane with the trajectory of an ion penetrating the surface at O. We 
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choose the X axis to lie in that plane. Finally, Y is the remaining direction 
which completes a right-handed reference frame, see Fig. [I]. 

(ii) Next we relate the quantities measured in coordinates of the local 
frame to coordinates in the laboratory frame (x,y,h). The latter is defined 
by the experimental configuration. That is, h is the direction normal to the 
uneroded flat surface. The ion trajectories together with the h axis define a 
plane, which is taken to be the x — h plane. And finally the y axis completes 
a right-handed reference frame, see Fig. |l|. However, <p, which is the angle 
between the ion trajectory and the local normal to the surface, changes from 
point to point along the surface, and is a function of the local values of the 
slopes at O (as seen in the laboratory frame) , as well as of the fixed angle 9 
subtended by the ion trajectories and the normal to the uneroded surface (the 
h direction in Fig. [l]). 

(iii) In the absence of overhangs the surface can be described by a single 
valued height function h(x,y,t), measured from an initial flat configuration 
which is taken to lie in the (x,y) plane. The ion beam is parallel to the x-h 
plane forming an angle < 9 < 7r/2 with the z axis. To obtain the equation 
of motion for the surface profile function h{x, y, t), we will have to project the 
normal component of the velocity of erosion onto the global h axis. We get 

dh{X d ^ t] — —v((p, Rx-, Ry) VTTW, (7) 

where ip is the angle of the beam direction with the local normal to the surface 
at h(x, y, t) and Rx,y the values of the local radii of curvature at (x, y, h). Now 
ip is a function of the angle of incidence 9 and the values of the local slopes 
d x h and d y h, and can be expanded in powers of the latter. We will assume 
that the surface varies smoothly enough so that products of derivatives of h 
can be neglected for third or higher orders. 

At this stage additional relevant physical processes must be taken into 
account to describe the evolution of the surface. First, the bombarding ions 
reach the surface at random positions and times. We account for the stochastic 
arrival of ions by adding to (^) a Gaussian white noise rj(x,y,t) with zero 
mean and variance proportional to the flux J. Second, at finite temperature 
atoms diffuse on the surfaceflO. To include this surface self-diffusion we allow 
for a term — DtV 2 (V 2 /i)E3'c3, where Dt is a temperature dependent positive 
coefficient. Expanding (|7|) and adding the noise and the surface-diffusion terms 
we obtain the equation of motionEH 



dh dh d 2 h d 2 h X T ( dh\ z X y f dh 
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From ([?]) we can compute the expressions for the coefficients appearing in (||) 
in terms of the physical parameters characterizing the sputtering process. To 
simplify the discussion we restrict ourselves to the symmetric case a = fi. The 
general case is discussed inEZl. If we write F = (e,Jp/\f2Ti) exp(— a 2 /2 — a 2 s 2 ), 
s = sin#, c = cos9 and a a = a/a, we find for the coefficients in (^) 

,2„2 i\ 



v = 


F F 
— c , 7 = — s{( 
a a 


X x — 


F ( 2 to 2 5 

—c\a a (3s —c 
a 


Xy = 


p 

c{a 2 c 2 }, 

a 


v x = 


^a a {2s 2 -c 2 - 


Vy = 


F 2 


Dl = 


2Aa a 



4 2 21 

a„s c 



2 2 2"! 
a a s C ) , 



(9) 



5 Analysis of the obtained growth equations 

Consistent with the direction of the bombarding beam and the choice of coor- 
dinates, the terms in (^) are symmetric under y — > — y but not under x — > —x, 
while for 9 — ► CL-we get 7 = £ x = ^ y = 0, A^ = X y and i/ x — v y . The equation 
studied in Ref . c3 corresponds to the deterministic linear version of (|J), i. c. 
\ x = \ y =rj = 0. 

If v x and v y are positive, the surface diffusion term is expected to con- 
tribute negligibly to the relevant surface relaxation mechanism when we probe 
the system at increasingly large length scales. Scaling properties are then de- 
scribed by the anisotropic KPZ equation (AKPZ), which predicts two possible 
behaviors depending on the relative signs of the coefficients X x and X y on. If 
A x Aj, > 0, then a = 0.38 and z = 1.6, the surface width W{L,t) increases al- 
gebraically—being characterized by the exponents of the KPZ equation in 2+1 
dimensionsE9. For A^Ay < 0, the nonlinear terms X x and X y become irrelevant, 
and the width grows only logarithmically, i.e. a = 0. 
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In our case v x can change sign as 8 and a a are varied, while v y is always 
negative. The negative v causes an instability, whose origin is the faster erosion 
for the bottom of a trough than for the peak of a crest, as predicted by (H)l3E] 
(see also Fig. 3 of Ref. c3) . An instability due to a negative surface tension 
is also known to take place in the Kuramoto-Sivashinsky (KS) equation l3, 
which is the noiseless and isotropic version of (g). It has been argued Jbr the 
KS equation that in 1+1 dimensions v renormalizes to a positive valueEil, and 
the large length scale behavior is described by the KPZ equation. In 2+1 
dimensions it is not completely settled whether the large distance behaviors of 
KS and KPZ fall in .the same universality class, different approaches leading 
to conflicting resultsl23. 

In contrast to the KS equation, Eq. (^) is anisotropic, and explicitly con- 
tains a noise term. The competition between surface tension and surface dif- 
fusion generates a characteristic length scale in the system, £ c = tJDt/\v\, 
where v is the largest in absolute value of the negative surface tension coeffi- 
cients. Below we discuss a possible scenario for the scaling behavior predicted 
by (|J) based primarily on the results available in the literature for some of 
its limits. The complete scaling picture should be provided by either a DRG 
analysis capable of coping with the linear instabilities present in the system, 
or a numerical integration of (^). 

The scaling behavior depends on the relative signs of u x , v yi X x and Aj,El. 
The variations of these coefficients as functions of a a and 8 lead to the phase 
diagram shown in Fig. ||. 

Regions I and IL=r For small 8 both v x and v y are negative. As discussed by 
Bradley and HarpertJ and experimentally studied by Chason et al.a, a periodic 
structure dominates the surface morphology, with ripples oriented along the 
direction (x or y) which presents the largest absolute value for its surface 
tension coefficient. The observed wavelength of the ripples is A c = 2tt\/2£ c . 

The large length scale behavior £ ^> £ c is expected to be different. Now 
both nonlinearities and the noise may become relevant. The scaling properties 
of the surface morphologies predicted by (||) are unknown. A possible scenario 
is that the u's renormalize to positive values, as they do for the KS equation 
in 1 + 1 dimensions, and the large scale scaling properties of the system are 
described by the AKPZ equation. Then one would observe algebraic scaling 
in region I, where both nonlinearities have the same (negative) sign, whereas 
scaling would become logarithmic through an AKPZ-like mechanism in region 
II, where and \ y have opposite signs. Actually, the asymptotic KPZ scal- 
ing has been recently shown to. occur along the 8 axis of Fig. 2 through a 
renormalization group analysisE-3. 

Region III — This region is characterized by a positive v x and a negative 
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Figure 2: Phase diagram for the isotropic case a ■■ 



1. Region I: v x < 0, v y < 0, A^ < 0, 



X y < 0; Region II: u x < 0, v y < 0, A^ > 0, X y < 0; Region III: u x > 0, v y < 0, A^ > 0, 
X y < 0. Here a is measured in arbitrary units and 6 is measured in degrees. 



v y . Now the periodic structure associated with the instability is directed along 
the y direction and is the dominant morphology at scales t ~ t c . Again, such 
an anisotropic and linearly unstable equation is unexplored in the context of 
growth equations. Assuming that v y renormalizes to a positive value, and that 
the AKPZ mechanism operates, one would expect logarithmic scaling in region 
since the nonlinear terms have opposite signs. 

Even though several aspects of the scaling behavior predicted by (8) and 
(^) remain to be clarified, we believe that these equations contain-the relevant 
ingredients for understanding roughening by ion bombardment To sum- 
marize, at short length scales the morphology consists of a periodic structure 
oriented along the direction determined by the largest in absolute value of the 
negative surface tension coefficients Modifying the values of a a or 9 changes 
the orientation of the ripplesBa. At large length scales we expect two different 
scaling regimes. One is characterized by the KPZ exponents, which might be 
observed in region / in Fig. ^|. Indeed, the values of the exponents reported by 
Eklund et al. li-3 are consistent within the experimental errors with the KPZ ex- 
ponents in 2+1 dimensions. The other regions {II and III) are characterized 
by logarithmic scaling [a = 0), which has not been observed experimentally so 
far. Moreover, by tuning the values of 9 and/or a a one may induce transitions 
among the different scaling behaviors. For example, fixing a a and increasing 
the value of 9 would lead from KPZ scaling (region J) to logarithmic scaling 
(//, 777) for large enough angles. 

Recent results by Rost and Krug on the two dimensional anisotropic KS 
equation indicate that the scaling regimes II and III in the noiseless limit of 
our model is dominated by exponentially growing solutions of the KS equation 
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. In those regions the ripple structure is oriented along a direction which is 
neither x nor y. Further numerical simulations are needed to understand the 
effect of the noise on the stability of the exponential solutions. Insight into 
the expected morphologies is obtained from numerical simulation of discrete 
models, correctly capturing the basic mechanisms taking place during sputter- 
ing. Recent simulations on discrete models indicate that the noisy KS equation 
indeed describes the dynamics of the sputtering generated roughening l3. 

6 Ion-Induced Surface Diffusion in Ion Sputtering 

In the absence of ion bombardment surface diffusion is thermally activated, and 
characterized by the diffusion constant, Dt — Dq exp [—Ed/kgT], such that 
the evolution of the surface ifiight, h(x,y,t) is described by the continuum 
equation dh/dt = — £>rV 4 /iE3. Here is the activation energy for surface 
diffusion of the adatoms and T is the substrate temperature. However, numer- 
ous experiments regarding the effect of ion bombardment on island formation, 
surface migration, surface smoothing and ripple formation have provided evi- 
dence that ion bombardment is accompanied by an increase in surface diffusion 
^C3>He1Ic3'c2L In particular, it has been demonstrated-that ion-induced sur- 
face diffusion caii decrease the epitaxial temperature E3, enhance nucleation 
during growth ej, modify the surface morphology, or induce the existence of 
a temperature independent ion-induced surface diffusion constant, as in the 
experiments of MacLaren et al. referred to above. 

Although the effect of the ions on surface diffusion is well documented 
experimentally, there is no theory that would quantify it. Eq. (|^) provides 
analytically the ion- induced diffusion constant, D 1 , and its dependence on the 
ion energy, flux, angle of incidence, and penetration depth. Consistent with 
symmetry considerations for 8 = we obtain D\. = Dy. However, for 8^0 wc 
find that D' x ^ D*, i.e. the ion-induced surface diffusion is anisotropic. More- 
over, its sign also depends on the experimental parameters. Their properties 
can be summarized as follows: (a) Independent of the angle of incidence Dy 
is positive, and decreases with 9, while the sign of the D\. depends on both 
8 and a a as shown in Fig. |§|. Thus, while for 8 = the ion bombardment 
enhances the surface diffusion (D^ > 0), for large 8 it can suppress diffusion; 
(b) The fact that D\. can be negative indicates that any simple theory con- 
necting the magnitude of the ion-induced diffusion to the energy transferred 
by the ions to the surface is incomplete, since it can predict only a positive 
D 1 . In fact, D 1 is the result of a complex interplay between the local sur- 
face topography and the energy transferred to the surface; (c) The diffusion 
constants are proportional to the flux J, in agreement with the detailed exper- 
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Figure 3: Ion-induced diffusion constant, D\. and Dy (inset) as a function of the angle of 
incidence 8. In both figures the curves correspond to a a = 1.5 (dashed line), a a = 2.0 
(dotted line) and a a = 2.5 (continuous line). 



imental study of Cavaille and Dreschnerli3; (d) It is a standard experimental 
practice to report the magnitude of the ion-enhanced diffusion using an effec- 
tive temperature T e ^ at which the substrate needs to be heated to obtain the 
same mobility as with ion bombardment We can calculate T e ^ using the 
relation D 1 + Dq exp(—E a /kBT) = Dq ex-p(—E a /kBT e ''), that has two impor- 
tant consequences. First, the anisotropic diffusion constant translates into an 
anisotropic T e ^ , i.e. we have ^ TyM . The experimental methods used 
to estimate T e ^ could not distinguish and T*ff c3c3. However, current 
observational methods should be able to detect the difference between the two 
directions. Second, while it is generally believed that ion bombardment can 
only raise the effective temperature since it transfers energy to the surface, the 
negative D 1 indicates that along the x direction one could have T e ff (e) 
Finally, the results (||) are based on Sigmund's theory of sputtering E3, that 
describes sputtering in the linear cascade regime. The energy range when this 
approach is applicable lies between 0.5 keV and 1 MeV, the precise lower and 
upper limits being material dependent. Thus, we do not expect (Q) to apply 
to low energy (few eV) ion-enhanced epitaxy. 

Quantitative comparison with experiments — At nonzero temperature the 
total diffusion constant is given by D = D 1 + Dt- As T decreases there 
is a critical temperature, T c , at which D 1 = Dt, so that for T < T c the 
diffusion is dominated by its ion-induced component, which is independent of 
temperature, in agreement with the experimental results of MacLaren et al. 
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0. Unfortunately, for most materials the quantities entering in F, a a and 
Do are either unknown, or only their order of magnitude can be estimated. 
However, we can express T c in terms of measurable quantities independent of 
these constants 

I 1 

c l-(2T k B /E a )ln(e Ion /e To ) K ' 

where £t is the experimentally measured ripple wavelength at any temperature 
To > T c ; lion is the ripple wavelength in the low temperature regime, T < T c , 
where ion induced diffusion dominates, and therefore £j on is independent of 
T; E a is provided by the slope of ln(^) versus 1/T in the high temperature 
regime (T >> T c ). Consequently, all quantities in ( p^| ) can be obtained from 
a plot of the ripple wavelength as a function of temperature, so (10) gives T c 

S terms of measurable quantities. Such a plot is provided by MacLaren et al. 
, leading to E a = 0.51eV, £ Ion = 0.8/zm. Using £ To = 2/zm for T = 368K, 
we obtain T c = 57° C, which is irugood agreement with the experiments, that 
provide T c between 45 and 60° CB. 

Morphological phase diagram — The detailed morphological phase diagram 
is rather complex if the diffusion is not thermally activated, but ion-induced. 
At low temperatures, when Dt is negligible, the ripple wavelengths are i\. = 



2ix^/D x j\v x \ and £ y = 2irJ Dy/\v y \. In the following we discuss the expected 

surface morphologies in function of the experimental parameters and a a , 
based on the phase diagram shown in Fig. [I| 

Region I — The surface tensions, v x and v y , are negative, while D x and D y 
are positive, consequently we have a superimposed ripple structure along the 
x and the y directions. The ripple wavelength observed experimentally is the 
smallest of the two, and since l\. > £ y the ripple wave vector is oriented along 
the y direction. The lower boundary of this region separating it from Region 
II is given by the solution of the i\. = £ y equation. 

Region II — Here the ripple wave vector is oriented along the x direction, 
since 4 < £\. This region is bounded below by the D\. = line. At large 
length spQ^sJpj^.nH II one expects kinetic roughening described by the KPZ 

Region III — In this region D x is negative, while the sign of all other 
coefficients are as in I and II. Since both the surface tension and the surface 
diffusion are destabilizing along x, every mode is unstable ancLone expects 
that the KPZ nonlinearity cannot turn on the KS stabilizationEZlo, the system 
being unstable at large length scales as well, leading to exponential growth. 
The lower boundary of this region is given by the v x = line. 

Region IV — Here we have v x > 0, v y < 0, D x < and D y > 0, i.e. one 
expects the surface to be periodically modulated in the y direction. In the x 



equation! 
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Figure 4: Phase diagram for the isotropic case a = fx = 1 at T < T c . Region I: v x < 0, 
u y < 0, D x > 0, D' y > and t x > t y ; Region II: v x < 0, u y < 0, D x > 0, D y > 0, and 
i x <l v \ Region III: ia, < 0, v y < 0, < and D y > 0; Region IV: v x > 0, v y < 0, < 
and > 0. Note that the phase diagram is independent of the precise values of J and p, 
while the e dependence is contained in <v. 

direction we have an interesting reversal of the instability: the short length 
scale instability generated by the negative D T X is stabilized by the positive 
surface tension v x . Thus there is no ripple structure along the x direction. 
Regarding the large length scale behavior, along the x direction the surface 
diffusion term is irrelevant compared to the surface tension, thus one expects 
KPZ scaling. However, along the y direction the KS mechanism is expected 
to act, renormalizing the negative v y to positive values for length scales larger 
than £y, leading to a large wavelength KPZ behavior. 

If thermal and ion-induced diffusion coexist, the ripple wavelengths are 



given by 4 = 2tt[(£> t + D* x ) /K|] 1/2 and t y = 2tt[(D t + D 1 )/]^]} 1 / 2 . The 



phase diagram for intermediate temperatures can be calculated using the total 
D. In particular for high T, when >> OL and » D, 1 , the phase 



diagram converges to the one obtained in Ref .113 the ripple orientation being 
controlled by the v x — v y line (dotted line in Fig. ||). Thus with increasing Dt 
the phase boundary between the regions I and II converges to the v x = v y line 
and the D x = boundary separating the regions II and III shifts downwards, 
eventually disappearing. However, in the intermediate regions new phases with 
coarsening ripple domainsES appear as the D x = line crosses the v x = line. 

While we limited our discussion to the effect of the ion-induced diffusion 
on the surface morphology, the results (|h can be used to investigate other 
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phenomena as well, such as island nucleation. The experimental verification of 
the above results would constitute an important step to elucidate the mecha- 
nism responsible for ion-induced diffusion, with potential applications to ion- 
enhanced epitaxy as well. 

7 Atomistic models 

The methods employed for the modeling of growth phenomena at the atomic 
level range from first-principle calculations, to molecular dynamics, and Monte 
Carlo (MC) simulations ej. However, currently only MC methods can reach 
the long time scales and fairly large length scales needed to observe ripple for- 
mation and roughening. Moreover, the general features of ripple formation and 
roughening are rather generic, suggesting the existence of a material indepen- 
dent robust mechanism governing them. Thus simple models that incorporate 
the basic physics of the system, i.e. ion-bombardment and surface diffusion, 
should be successful in capturing the observed behavior, see e.g. Ref. EH. We 
have developed atomistic MC models of erosion that include ion-bombardment, 
ion-induced atom removal and activated surface diffusion. To model surface 
diffusion we use the methods developed for MBE, taking the diffusion rate of 
an atom proportional to exp[— Ea/ksT], where E a is the activation energy for 
diffusion and T is temperature!!!. 

To model the ion-bombardment we assume that the ion beam has a con- 
stant flux, but the time and the position when and where an ion strikes the 
surface is random. The ion penetrates the bulk reaching a penetrating depth 
a and releases its kinetic energy (see inset of Fig. 1). The energy reaching 
the surface atoms, Ei on , is calculated using the energy distribution (5). If 
Ei on is larger than the desorption energy Ed = nEs + Eq + Ed, where n is 
the number of nearest neighbors, then it will be sputtered. If Ei on is smaller 
than Ed, then it will contribute either to surface diffusion with probability 
exp[— (jiEb + Eq — Ei on )/kBT] or to breaking the bonds (sputtering) with 
probability exp[—(Ed — Ei on ) /fesT]. In the simulation we used e=1000 eV, 
<7=/i=a=10.9 lattice constants, Eq = 0.4 eV, Eb = 0.1 cV, and Ed=l-0 eV. 

We used the structure factor, S(k) =< h(k,t)h(—k,t) >, where h(k,t) is 
the Fourier transform of h{x, t), to estimate the ripple wavelength I. Typically 
S(k) develops a sharp maximum, which allows us to estimate £. However, at 
large erosion times the structure factor reflects the roughening of the substrate 
as well, thus the maximum will slowly disappear. Similarly, the maximum is 
more visible at low temperatures than at large T. 

Fig. P shows the ripple wavelength as a function of temperature for 9 = 0. 
At large temperatures I follows an Arrhenius law, however, at low T it con- 
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Figure 5: Ripple wavelength as a function of the inverse temperature obtained in the nu- 
merical simulations of a one dimensional sputtering model. The continuous line represents 

an analytical fit (see the text). 



verges to a constant value. These results are reminiscent of the experimental 
data of Maclaren et al. Ej, providing direc t numerica l evidence of ion enhanced 
surface diffusion. Indeed, using I = 2tt\/ Dt + D 1 , where Dt follows an Ar- 
rhenius law Dt = Doex-p(—E a /kBT), we can obtain an excellent fit to the 
data in Fig. ||. Furthermore, the fit also allows us to determine the diffu- 
sion coefficients in our simulation, providing Do = 7000, E a = 0.08eV, and 
D 1 = 1543. There are two observations we have to make analyzing these re- 
sults. First, E a is much lower than the activation energy we used as an input 
to the simulations. However, this is not in contradiction: the effective activa- 
tion energy felt by the diffusing atoms is lowered by the energy provided by 
the ions. This effect lowers E a . Indeed, we measured the average activation 
energy E%?f = E a — Ei for the atoms, the result agreeing in order of mag- 
nitude with the obtained from the fit in Fig. ||. Second, D is smaller 
than the experimentally expected values, that are of order 10 xx , but agrees 
with the smaller diffusive activity we can obtain in the numerical simulations 
(due to running time limitations). It is easy to see that a larger Do would 
lead to a more extended high-temperature region, such as the one observed by 
MacLaren et al. E3. 
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